The package OpenPNM allows for the simulation of many transport phenomena in porous media such as Stokes flow, Fickian diffusion, advection-diffusion, transport of charged species, etc. Transient and steady-state simulations are both supported. An example of a transient Fickian diffusion simulation through a Cubic pore network is shown here.
First, OpenPNM is imported.
In [1]:
import numpy as np
import openpnm as op
np.random.seed(10)
%matplotlib inline
np.set_printoptions(precision=5)
In [2]:
ws = op.Workspace()
ws.settings["loglevel"] = 40
proj = ws.new_project()
In [3]:
net = op.network.Cubic(shape=[29, 13, 1], spacing=1e-5, project=proj)
Here, a geometry, corresponding to the created network, is created. The geometry contains information about the size of pores and throats in the network such as length and diameter, etc. OpenPNM has many prebuilt geometries that represent the microstructure of different materials such as Toray090 carbon papers, sand stone, electrospun fibers, etc. In this example, a simple geometry known as StickAndBall that assigns random diameter values to pores throats, with certain constraints, is used.
In [4]:
geo = op.geometry.StickAndBall(network=net, pores=net.Ps, throats=net.Ts)
In [5]:
phase = op.phases.Water(network=net)
Next, a physics object is defined. The physics object stores information about the different physical models used in the simulation and is assigned to specific network, geometry and phase objects. This ensures that the different physical models will only have access to information about the network, geometry and phase objects to which they are assigned. In fact, models (such as Stokes flow or Fickian diffusion) require information about the network (such as the connectivity between pores), the geometry (such as the pores and throats diameters), and the phase (such as the diffusivity coefficient).
In [6]:
phys = op.physics.GenericPhysics(network=net, phase=phase, geometry=geo)
The diffusivity coefficient of the considered chemical species in water is also defined.
In [7]:
phase['pore.diffusivity'] = 2e-09
In [8]:
mod = op.models.physics.diffusive_conductance.ordinary_diffusion
phys.add_model(propname='throat.diffusive_conductance', model=mod, regen_mode='normal')
In [9]:
fd = op.algorithms.TransientFickianDiffusion(network=net, phase=phase)
In [10]:
fd.set_value_BC(pores=net.pores('front'), values=0.5)
fd.set_value_BC(pores=net.pores('back'), values=0.2)
In [11]:
fd.set_IC(0.2)
Note that both set_value_BC and set_IC also accept as input, in addition to a single scalar value, an ndarray.
The settings of the transient algorithm are updated here. This step is optional as default settings are predefined. It is, however, important to update these settings on each new simulation as the time-scale of different phenomena in different problems may strongly differ.
Here, the time discretization scheme is set to cranknicolson, which is second-order accurate in time. The two other options supported in OpenPNM are the implicit scheme (only first order accurate but faster than the cranknicolson) and the steady which simply corresponds to a steady-state simulation.
Other parameters are also set; the final time step t_final, the output time stepping t_output, the computational time step t_step, and the tolerance to be achieved before reaching steady-state t_tolerance.
In [12]:
fd.setup(t_scheme='cranknicolson', t_final=100, t_output=5, t_step=1, t_tolerance=1e-12)
Note that the output time stepping t_output may a scalar, ND-array, or list. For a scalar, it is considered as an output interval. If t_output > t_final, no transient data is stored. If t_output is not a multiple of t_step, t_output will be approximated. When t_output is a list or ND-array, transient solutions corresponding to this list or array will be stored. Finally, initial, final and steady-state (if reached) solutions are always stored.
In [13]:
print(fd.settings)
Note that the quantity corresponds to the quantity solved for.
In [14]:
fd.run()
In [15]:
print(fd)
Note that the solutions at every exported time step contain the @ character followed by the time value. Here the solution is exported after each $5s$ in addition to the final time step which is not a multiple of $5$ in this example.
To print the solution at $t=10s$
In [16]:
fd['pore.concentration@10']
Out[16]:
The solution is here stored in the phase before export.
In [17]:
phase.update(fd.results())
Export the results into an xdmf file to be able to play an animation of the time dependent concentration on Paraview.
In [18]:
proj.export_data(phases=[phase], filename='./results/out', filetype='xdmf')
In [19]:
#NBVAL_IGNORE_OUTPUT
import matplotlib.pyplot as plt
c = fd['pore.concentration'].reshape((net._shape))
fig, ax = plt.subplots(figsize=(6, 6))
plt.imshow(c[:,:,0])
plt.title('Concentration (mol/m$^3$)')
plt.colorbar();